;**************************************************************
; Calculate the "joint PDF [%]" of arrays "x" and "y"
;**************************************************************
undef("pdfxy_ylog")
function pdfxy_ylog(x:numeric, y:numeric, nbinx[1]:integer, nbiny[1]:integer, opt:logical)
local nGood_x, nbinsx, xMin, xMax, mnmxint, binx, binxBound, binxCenter, nbx \
     ,nGood_y, nbinsy, yMin, yMax, ySpace , biny, binyBound, binyCenter, nby \
     ,outside, mnmxint, pdf2, nTot, dimx, dimy, rankx, ranky, x1d, y1d, iy, warnFlag\
     ,xMIN, xMAX, yMIN, yMAX, binxBoundMin, binxBoundMax, binyBoundMin, binyBoundMax\
     ,nMax, nUse, nxLoOut, nxHiOut, nyLoOut, nyHiOut, popt, fmsg
begin

  fmsg = default_fillvalue("float")

  if (nbinx.le.2) then
      nbinsx = 25             ; default
  else
      nbinsx = nbinx          ; user sprecified
  end if

  if (nbiny.le.2) then
      nbinsy = 25             ; default
  else
      nbinsy = nbiny          ; user specified
  end if
                              ; error check
  warnFlag = True             ; default ... print warning messages
  if (opt .and. isatt(opt,"WarningMsg") .and. .not.opt@WarningMsg) then
      warnFlag = False
  end if

  if (warnFlag) then
      nGood_x = num(.not.ismissing(x))
      if (nGood_x.lt.3) then
          print("pdfxy: nGood(x)="+nGood_x+" : Need more non-missing points")
      end if
    
      nGood_y = num(.not.ismissing(y))
      if (nGood_y.lt.3) then
          print("pdfxy: nGood(y)="+nGood_y+" : Need more non-missing points")
      end if
    
      if (nGood_x.lt.3 .or. nGood_y.lt.3) then
          print("pdfxy: exit error: not enough points")
          pdf2 = new( (/nbinsy,nbinsx/), "double", 1d20)
          return( pdf2 )
      end if
  end if              ; warnFlag

  xMIN = min(x)*1d0   ; min/max for ENTIRE array             
  xMAX = max(x)*1d0   ; force "double" [*1d0]
  yMIN = min(y)*1d0   
  yMAX = max(y)*1d0  
                                             ; max number possible
  nMax = num(x.ge.xMIN .and. x.le.xMAX .and. \
             y.ge.yMIN .and. y.le.yMAX) 

  if (nMax.lt.3) then                        
      print("pdfxy: exit error: nMax.lt.3")
      pdf2 = new( (/nbinsy,nbinsx/), "double", 1d20)
      return( pdf2 )
  end if

  xMin = 0.0d         ; prototype as "double"    
  xMax = 0.0d
  yMin = 0.0d     
  yMax = 0.0d
                      ; User may want different min/max used
  if (opt .and. isatt(opt,"binx_min")) then  
      xMin = opt@binx_min                   ; user
  else
      xMin = xMIN                           ; calculated
  end if

  if (opt .and. isatt(opt,"binx_max")) then  
      xMax = opt@binx_max                   ; user
  else
      xMax = xMAX                                         
  end if

  if (opt .and. isatt(opt,"biny_min")) then  
      yMin = opt@biny_min                   ; user
  else
      yMin = yMIN
  end if

  if (opt .and. isatt(opt,"biny_max")) then  
      yMax = opt@biny_max                   ; user
  else
      yMax = yMAX  
  end if
                                            
  if (opt.and.isatt(opt,"binx_nice")) then  ; nice xMin, xMax
      outside = False                       
      if (isatt(opt,"binx_nice_outside")) then
          outside = opt@binx_nice_outside
      end if
      mnmxint = nice_mnmxintvl( xMIN, xMAX, nbinsx, outside)
      xMin    = mnmxint(0)
      xMax    = mnmxint(1)
      xSpace  = mnmxint(2)
    ;;nbinsx  = round( (xMax-xMin)/xSpace , 3) 
      delete(mnmxint)
  end if

  if (opt .and. isatt(opt,"biny_nice")) then ; nice yMin, yMax
      outside = False
      if (isatt(opt,"biny_nice_outside")) then
          outside = opt@biny_nice_outside
      end if
      mnmxint = nice_mnmxintvl( yMIN, yMAX, nbinsy, outside)
      yMin    = mnmxint(0)
      yMax    = mnmxint(1)
      ySpace  = mnmxint(2)
    ;;nbinsy  = round( (yMax-yMin)/ySpace , 3) 
      delete(mnmxint)
  end if

;;dbinx       = (xMax-xMin)/(nbinsx-1)
  dbinx       = (xMax-xMin)/nbinsx
  binxBound   = xMin + ispan(0,nbinsx,1)*dbinx
  binxBound(nbinsx) = xMax                       ; avoid roundoff
  binxCenter  = (binxBound(0:nbinsx-1) + binxBound(1:nbinsx))*0.5d 

;;dbiny       = (yMax-yMin)/(nbinsy-1)
  dbiny       = (yMax-yMin)/nbinsy
;;kaizhang 
  binyBound   = exp(fspan(log(yMin),log(yMax),nbinsy+1)) ;;yMin + ispan(0,nbinsy,1)*dbiny
  binyBound(nbinsy) = yMax                       ; avoid roundoff
  binyCenter  = (binyBound(0:nbinsy-1) + binyBound(1:nbinsy))*0.5d 

  binxBoundMin= binxBound(0)          ; convenience
  binxBoundMax= binxBound(nbinsx)   
  binyBoundMin= binyBound(0)   
  binyBoundMax= binyBound(nbinsy)   

  pdf2        = new( (/nbinsy,nbinsx/), "double", 1d20)
  pdf2        = 0.0d                  ; initialize
  
  dimx        = dimsizes(x)
  dimy        = dimsizes(y)
  rankx       = dimsizes(dimx)
  ranky       = dimsizes(dimy)

  popt          = True
  popt@fraction = False     ; False means to return %
  pdf2 = pdfxy_bin(x,y, binxBound,binyBound, popt)
                          ; all NCL [original code]
    ;;x1d  = ndtooned(x)  ; requires more memory
    ;;y1d  = ndtooned(y)

    ;;do nby=0,nbinsy-1
    ;;   iy = ind(y1d.ge.binyBound(nby) .and. y1d.lt.binyBound(nby+1))
    ;;   if (.not.ismissing(iy(0))) then
    ;;       do nbx=0,nbinsx-1
    ;;          pdf2(nby,nbx) = num(x1d(iy).ge.binxBound(nbx) .and. \
    ;;                              x1d(iy).lt.binxBound(nbx+1)     )
    ;;       end do
    ;;   end if
    ;;   delete(iy)      ; size may change
    ;;end do

  if (popt .and. isatt(popt,"fraction") .and. popt@fraction) then
      pdf2 = 1d2*pdf2                ; frequency [%]
  end if

  pdf2!0   = "y"         ; arbitrary name
  pdf2!1   = "x"
  pdf2&x   =  binxCenter 
  pdf2&y   =  binyCenter
                                             ; actual number used
  nUse     = num(x.ge.binxBoundMin .and. x.le.binxBoundMax .and. \
                 y.ge.binyBoundMin .and. y.le.binyBoundMax)

  if (nMax.eq.nUse) then
      nxLoOut  = 0
      nxHiOut  = 0
      nyLoOut  = 0
      nyHiOut  = 0
  else
      nxLoOut  = num(x.lt.binxBoundMin)
      nxHiOut  = num(x.ge.binxBoundMax)
      nyLoOut  = num(y.lt.binyBoundMin)
      nyHiOut  = num(y.ge.binyBoundMax)
  end if

  pdf2@nMax            = nMax
  pdf2@nUse            = nUse
  pdf2@nOut            = nxLoOut + nxHiOut + nyLoOut + nyHiOut
  if (nMax.gt.0) then
      pdf2@pcUse       = (1e2*nUse)/nMax         ; percent [pc]
      pdf2@pcOut       = (1e2*pdf2@nOut)/nMax
  else
      pdf2@pcUse       = fmsg
      pdf2@pcOut       = fmsg
  end if

  pdf2@xMIN            = xMIN 
  pdf2@xMAX            = xMAX 
  pdf2@binx_center     = binxCenter
  pdf2@binx_bounds     = binxBound
  pdf2@binx_bound_min  = binxBoundMin
  pdf2@binx_bound_max  = binxBoundMax
  pdf2@binx_spacing    = dbinx   ; binxBound(2)-binxBound(1)
  pdf2@nbinsx          = nbinsx
 ;pdf2@nxLoOut         = nxLoOut
 ;pdf2@nxHiOut         = nxHiOut
  
  if (nMax.gt.0) then
      pdf2@pcxLoOut    = (1e2*nxLoOut)/nMax         ; %
      pdf2@pcxHiOut    = (1e2*nxHiOut)/nMax
  else
      pdf2@pcxLoOut    = fmsg
      pdf2@pcxHiOut    = fmsg
  end if

  pdf2@yMIN            = yMIN 
  pdf2@yMAX            = yMAX 
  pdf2@biny_center     = binyCenter
  pdf2@biny_bounds     = binyBound
  pdf2@biny_bound_min  = binyBoundMin
  pdf2@biny_bound_max  = binyBoundMax
  pdf2@biny_spacing    = dbiny   ; binyBound(2)-binyBound(1)
  pdf2@nbinsy          = nbinsy
 ;pdf2@nyLoOut         = nyLoOut
 ;pdf2@nyHiOut         = nyHiOut
  if (nMax.gt.0) then
      pdf2@pcyLoOut        = (1e2*nyLoOut)/nMax         ; %
      pdf2@pcyHiOut        = (1e2*nyHiOut)/nMax
  else
      pdf2@pcyLoOut        = fmsg
      pdf2@pcyHiOut        = fmsg
  end if

 ;pdf2@long_name      = "Joint PDF"
 ;if (isatt(x,"long_name") .and. isatt(y,"long_name")) then
 ;    pdf2@long_name  = "Joint PDF: "+x@long_name+" | " \
 ;                                   +y@long_name
 ;end if
 ;pdf2@units          = "frequency"
  
  return( pdf2 )
end

